A polydisperse lattice-gas model 
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We describe a lattice-gas model suitable for studying the generic effects of polydispersity on 
liquid-vapor phase equilibria. Using Monte Carlo simulation methods tailored for the accurate de- 
termination of phase behaviour under conditions of fixed polydispersity, we trace the cloud and 
shadow curves for a particular Schulz distribution of the polydisperse attribute. Although poly- 
dispersity enters the model solely in terms of the strengths of the interparticle interactions, this is 
sufficient to induce the broad separation of cloud and shadow curves seen both in more realistic 
models and experiments. 



I. INTRODUCTION 

Lattice-based models of many-body systems have been 
in the vanguard of the drive for an improved understand- 
ing of equilibrium phase behaviour, almost since the in- 
ception of the field [lj . The utility of such models derives 
both from the simplifications to approximate analytical 
treatments that result from the spatial discretisation, and 
from their amenability to efficient numerical simulation. 
Indeed, even in the current era of teraflop computers, 
simulations of lattice models are still deployed routinely 
in fields as diverse as polymer physics and quantum chro- 
modynamics. 

The purpose of this paper is to introduce a new lattice- 
gas model suitable for simulation investigations of the 
phase behaviour of polydisperse fluids. These are com- 
plex fluids whose constituent particles exhibit continu- 
ous variation in terms of some physical attribute such as 
their size, shape or charge. Polydispersity is inherent to 
a host of natural and synthetic soft matter fluids includ- 
ing (inter alia) colloidal dispersions, polymer solutions, 
liquid crystals, as well as some biological fluids such as 
blood 2]. Understanding its generic effects on the phase 
behaviour of model systems is, therefore, not only a mat- 
ter of fundamental interest, but also one of considerable 
practical and commercial importance. 

The complexity bestowed on a fluid system by poly- 
dispersity complicates simulation studies of phase be- 
haviour. Although newly-developed methodologies [!, 0, 
H, @ allow the accurate study of phase coexistence in off- 
lattice models of polydisperse fluids, progress has never- 
theless been hindered (compared to studies of monodis- 
perse systems) by the fact that coexistence occurs over a 
region of the pressure-temperature phase diagram, rather 
than merely along a line [7], Q . The necessity of exploring 
a higher dimensional space of model parameters in order 
to fully characterize a given phase transition renders such 
investigations laborious. Hence progress in elucidating 
the generic effects of polydispersity on phase behaviour 
has been slower than be might be desired. 

Accordingly, there is a pressing need for a computa- 
tionally tractable model system, which, whilst not nec- 
essarily portraying realistically the microscopic features 
of any particular polydisperse fluid, does nevertheless 



capture principal features of the macroscopic phase be- 
haviour. The lattice-gas model presented below should 
go some way to meeting this need. Specifically, it should 
help answer unsolved questions such as how the coexis- 
tence and critical point properties depend on the nature 
of the particle interactions and the form and degree of the 
polydispersity. Additionally it should provide an efficient 
test bed for the development of fresh computational and 
analytical methodologies for studying polydisperse phase 
behaviour. 

The structure of our paper is as follows. In Sec. ILT1 we 
outline salient features of the phase behaviour of poly- 
disperse systems, before moving on to describe the poly- 
disperse lattice-gas (PLG) model in Sec. Mil We then 
report, in Sec. IIV1 a simulation study of key features 
of the liquid-vapor phase behaviour for one particular 
choice of the form of the polydispersity. Finally, Sec. [V] 
summarizes our findings and conclusions. 



II. BACKGROUND: POLYDISPERSE PHASE 
EQUILIBRIA 

In this section we present a brief overview of the princi- 
pal differences between phase coexistence in polydisperse 
and monodisperse systems, For a more detailed discus- 
sion, the interested reader is referred to a recent review 

The state of a polydisperse system (or any of its 
phases) is described by a density distribution pier), with 
p(o~)do~ the number density of particles whose polydis- 
perse attribute a lies in the range a . . . a+da. In the most 
commonly encountered experimental situation, the form 
of the overall or "parent" distribution, p^ \a), is fixed 
by the synthesis of the fluid, and only its scale can vary 
depending on the proportion of the sample volume occu- 
pied by solvent. One can then write p(°\a) = n^> f^°\a) 
where f^°'(o~) is the normalized parent shape function 
and = N/V the overall particle number density; as 
the latter is varied, p^(cr) traces out a "dilution line" 
in density distribution space. The phase diagram is then 
spanned by and the temperature T. 

The richness of the phase behaviour of polydisperse flu- 
ids stems from the occurrence of fractionation (To. ITU fl~2| : 
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at coexistence, particles of each a may partition them- 
selves unevenly between two or more "daughter" phases 
as long as-due to particle conservation-the overall den- 
sity distribution p^ (a) of the parent phase is main- 
tained. As a consequence, the conventional vapor-liquid 
binodal of a monodisperse system splits into a cloud curve 
marking the onset of coexistence, and a shadow curve giv- 
ing the density of the incipient phase; the critical point 
appears at the intersection of these curves rather than 
at the maximum of either Q. This splitting is seen in 
experiments on polydisperse fluids (see e.g. ref. [H[). 

The daughter distributions are related to the par- 
ent via a simple volumetric average: (1 — £) (a) + 
$,p^{cr) — p(°^(a) , where 1 — £ and £ are the fractional 
volumes of the two phases. In contrast to a monodisperse 
system where the densities of the coexisting phases are 
specified by the binodal and depend solely on tempera- 
ture, full specification of the coexistence properties of a 
polydisperse system requires not only knowledge of the 
cloud curve, but also the dependence of £, p^(cr) and 



2 )(cr) on and T @. 



III. POLYDISPERSE LATTICE-GAS MODEL 

Our polydisperse lattice-gas (PLG) model is defined 
within the grand canonical ensemble by the hamiltonian: 



H = - 



(aaTciWcjio-')-^ (1) 



Here a is the value (assumed scalar) of the polydisperse 
attribute, which controls the strength of interparticle in- 
teractions. We shall regard a as the label for a notional 
particle "species" , whose corresponding chemical poten- 
tial is p{o~). The exponent a can, in principle, take an 
arbitrary value but we shall restrict ourselves to the case 
a = 1 in the present study. Ci(o~) simply counts the 
number of particles of species a at site i, for which we 
impose a hard-core constraint such that X)cr c i( cr ) = 
or 1. The instantaneous density distribution follows as 
p{a) — L^ d '^2 i Ci(a), with d = 3 in the simulations re- 
ported below; i runs over the sites of a periodic lattice 
i — 1, L d , assumed simple cubic in this work. The sum 
in the first term on the right hand side of |T]) similarly 
runs over all pairs i,j of nearest neighbor sites, as well as 
over all combinations of a and a'. 

Put simply, each lattice site may be either vacant, or 
occupied by a single particle which carries a continuous 
species label a. Particles on nearest neighbor sites inter- 
act with a potential energy which is the negative of the 
product of their species labels. Thus the role of poly- 
dispersity in this model is to engender a spread of pos- 
sible interaction strengths between particles -a situation 
which contrasts with the single interaction str eng th char- 
acterizing the simple Ising lattice-gas model [14j . 

The model of Eq. Q] (which we have also briefly de- 
scribed in a different context elsewhere |5J|) has, to our 
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FIG. 1: The Schulz parent distribution for the case z — 50 
(cf. Eq. ©). 



knowledge, not been considered previously by other au- 
thors. We note that it is distinct from well-known lattice 
spin models such as the g-state Potts model in the limit 
of large q, and the XY model. The distinction arises 
both in terms of the form of the interactions, and the ab- 
sence of an isomorphic magnetic description, in contrast 
to the case of the simple (monodisperse) lattice-gas [3] . 



IV. SIMULATION STUDIES 

A. Parent distribution 

In this work we consider the case in which the species 
labels a are drawn stochastically from a (parental) dis- 
tribution of the Schulz form [151 ]: 



z + l 



z+l 



a exp 



z + l 



(2) 



with a mean species label a = 1. We have elected to 
study the case z = 50, corresponding to a moderate 
degree of polydispersity: the standard deviation of the 
species label is S = 1/ \Jz + 1 14% of the mean. The 
form of the distribution is shown in Fig.[TJ Although our 
motivation for employing the Schulz distribution is pri- 
marily ad-hoc, we note that it has been found to fairly 
accurately describe the polydispersity of some polymeric 
systems [16| . 

In the simulations reported below, the distribution 
f(°'(o~) was limited for computational convenience to the 
range 0.5 < a < 1.4, and renormalized appropriately. 
Imposition of lower and upper cutoffs, in the tails of the 
parent distribution, obviates the need to determine the 
associated portions of the chemical potential distribution 
p(a) conjugate to p(°\a) - a task that can be numeri- 
cally fraught when the parent density is very small. It 
should be noted, however, that depending on the form of 
the parent, the upper cutoff can have drastic effects on 
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the phase behaviour, even when it is located far into the 
tail of the distribution [171 . 



B. Methodology 

We have studied the phase behaviour of the PLG 
model, Eq. [TJ under conditions of fixed polydispersity 
via Monte Carlo simulations within the grand canonical 
ensemble (GCE). Our Metropolis algorithm invokes three 
types of lattice operations: particle deletions, particle in- 
sertions, and variation of a particle's species label a. The 
latter quantity is treated as a continuous variable in the 
permitted range 0.5 < a < 1.4; however, distributions 
defined on a, such as the instantaneous density p(a), 
and the chemical potential p(a), are represented as his- 
tograms formed by discretising this range into 100 bins. 
The simulation results presented below pertain to peri- 
odic cubic systems of linear dimension L = 10, 20, 30, 40 
lattice units. Further details concerning general aspects 
of the simulation of polydisperse fluids within the GCE, 
as well as the structure, storage and acquisition of data, 
have been presented elsewhere [|[ . 

The principal observable of interest is the fluctuat- 
ing form of the instantaneous density distribution p(a). 
From this we derive the distribution p(n) of the overall 
number density n = J dap(a). If one further identifies 
a with particle diameter - subject to the proviso that 
our model does not account for any diameter dependence 
of the hard core repulsion - we can consider, as an al- 
ternative to the density n, a notional "volume fraction" 
rj=(7r/6)Jdaa 3 p(a). 

The existence of phase coexistence at given chemical 
potentials is signalled by the presence of two distinct 
peaks in the probability distribution p(n). In order to 
obtain estimates of dilution line coexistence properties 
at some prescribed temperature, we employ an accurate 
approach recently proposed by ourselves [5|] . For a given 
choice of rv-°\ the method entails tuning the chemical 
potential distribution p,{a) together with a parameter £, 
such as to simultaneously satisfy both a gen eralized lever 
rule and an equal peak weight criterion [l8| for p{n): 



,(0)f(0) 



(l-eV (1 V)+£p (2 V) (3a) 
r = 1 (3b) 

Here the daughter density distributions p^(a) and 
p^ (a) are assigned by averaging only over configura- 
tions belonging to either peak of p(n) . The quantity r is 
the peak weight ratio: 



fn>n* P(n)dn 
/«<«• P^ dn 



(4) 



with n* a convenient threshold density intermediate be- 
tween vapor and liquid densities, which we take to be 
the location of the minimum in p(n). The tuning of p{<r) 



and £ necessary to simultaneously satisfy Eqs. (|3a[) and 
(|3b[) can be efficiently achieved by histogram extrapola- 
tion techniques pj| . given an initial estimate for p(a) 
which is provided by an iterative refinement technique, 
as described elsewhere 

The value of £ resulting from the application of the 
above procedure is the desired fractional volume of the 
liquid phase at the nominated value of . Cloud 
points are determined as the value of at which £ 
first reaches zero (vapor branch) or unity (liquid branch) , 
while shadow points are given by the density of the co- 
existing incipient daughter phase, which may be simply 
read off from the appropriate peak position in the cloud 
point form of p(n) and p(rf) . It should be pointed out 
that the finite-size corrections to estimates of coexistence 
properties obtained using the equal peak weight criterion 
for p{n) are exponentially small in the system size [5| . 

In order to obtain the phase behaviour of our model 
system, we scanned the dilution line for a selection of 
fixed temperatures. We started by setting T = T c , 
the critical temperature (known from a brief preliminary 
study ||), and tracked the locus of the dilution line in a 
stepwise fashion. This tracking procedure must be boot- 
strapped with knowledge of the form of fi(cr) at some 
initial point on the dilution line. A suitable estimate was 
obtained, for a point near the critical density, by means 
of the iterative refinement procedure [4], in conjunction 
with the equal peak weight criterion for p(n) discussed 
above. Simulation data accumulated at this near-critical 
state point was then extrapolated to a lower, but nearby 
density by means of histogram reweighting, thus pro- 
viding an estimate of the corresponding form of p(cr). 
The latter was employed in a new simulation, the re- 
sults of which were similarly extrapolated to a still lower 
value of n(°) . Iterating this procedure thus enabled the 
systematic tracing of the whole dilution line and the iden- 
tification of cloud points. Histogram extrapolation fur- 
ther permitted a determination of dilution line properties 
at adjacent temperatures, thereby facilitating a system- 
atic determination of the phase behaviour in the — T 
plane, including special loci such as the cloud curve and 
the critical isochore (n^ = n^). Im plem entation of 
multicanonical preweighting techniques |2(| at each co- 
existence state point ensured adequate sampling of the 
coexisting phases in cases where they are separated by a 
large interfacial free energy barrier. A fuller account of 
this latter procedure has been presented previously [2l| . 



C. Results 

The techniques outlined above were used to obtain ac- 
curate estimates of key features of the phase diagram of 
the PLG, namely the cloud and shadow curves, and the 
coexistence properties on the critical isochore. 

The liquid-vapor critical point has been located ap- 
proximately in a previous brief study of the PLG [5J, 

which found (ni°\T c ) — (0.521, 1.171) in reduced units. 
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This is to be compared with the critical parameters of 
the monodisperse (Ising) lattice-gas (0.5,1.127955) [22I ]. 
Thus the inclusion of polydispersity of the form (fl]) is 
seen to raise both the critical temperature and the critical 
density. Moreover we find that it splits the liquid-vapor 
binodal into well-separated cloud and shadow curves (cf. 
Fig. [2]) in a manner similar to that occurring in contin- 
uum fluid models for which polydispersity affects the in- 
terparticle interaction strength as well as its range [rj[l7l]. 
(When only the range is polydisperse, a scale symmetry 
forces the critical point to be near the top of cloud and 
shadow [H, [24], [2j| Hfjj].) As a consequence, the crit- 
ical point lies below the maximum of the cloud curve 
and phase coexistence can be observed even at T = T c , 
provided that < . Note that our procedure for 
cloud/shadow curve determination breaks down in the 
vicinity of the critical point, because the two peaks in 
p(n) overlap in a finite-sized system. This is the reason 
for absence of data near criticality in Fig. [21 
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curve well away from criticality. One sees that in keeping 
with previous findings for off-lattice models [f| [llj , there 
is considerable fractionation. Specifically, the number of 
particles in the liquid phase having a large value of a is 
strongly enhanced relative to the vapor phase. This en- 
hancement arises because particles having larger a inter- 
act more strongly than those of small a, thereby yielding 
a substantial free energy gain of the liquid (with its larger 
number of nearest neighbor particle contacts) over the 
vapor. One further sees that as a result of this enhance- 
ment, the form of the liquid daughter phase distribution 
is strongly affected by the presence of the cutoff. The 
issues surrounding such "cutoff effects" and their impli- 
cations for phase behaviour have recently been examined 
in detail in refs. ml 03 • 
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FIG. 3: Normalized daughter distributions f v (c) and fi(a) for 
the vapor (v) and liquid (I) phases as obtained for the L = 30 
system on the cloud curve at = 0.1275, T = 1.111, as 
described in the text. 

Finally, in this section, we have determined the co- 
existence curve corresponding to a parent density which 
is fixed at its critical value nf' = 0.521. This curve is 
presented in Fig. [4] One notes that its general shape 
is similar to the standard binodal that one finds in a 
monodisperse fluid. This is because the fractional vol- 
umes of the two phases are approximately (though not 
strictly) equal on this coexistence curve. By contrast, 
on the cloud curve the fractional volume of one phase 
is infinitesimal, leading to much stronger fractionation 
effects. 



V. CONCLUSIONS 



FIG. 2: Estimates of the cloud and shadow curves of the PLG 
model for the Schulz parent distribution of Eq. [2] obtained 
as described in the text. Data are shown for three system 
sizes, (a) Density-temperature plane, (b) Volume fraction- 
temperature plane. In each case the estimated critical point 
is marked by a star (*). 

Fig. [3] shows the normalized forms of the vapor and liq- 
uid daughter phase distributions for a point on the cloud 



In this paper, we have introduced a lattice-gas model 
for a polydisperse fluid and shown that it exhibits quali- 
tatively similar liquid-vapor phase behaviour to realistic 
continuum fluid models in which polydispersity affects 
both the size and the strength of the particle interactions 
d, [TtJ ■ In tests we find that the coexistence properties of 
the PLG can be determined with a computational effort 
that is smaller by a factor of 1 — 2 orders of magnitude 
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of the coexistence curve for which the 
parent density n {u> is fixed to its critical point value ie. 

Data is shown for three system sizes. The marked 
densities at each temperature derive from the peak positions 
in the corresponding form of p(n). Clear finite-size effects are 
visible in the vicinity of the critical point. 



compared to that for an off-lattice system having a simi- 
lar number of particles at equivalent state points such as 
criticality. This bodes well for the prospects of harnessing 
the PLG model to elucidate how coexistence properties 
are affected by factors such as the form of the parent 
distribution, the choice of the large-cr cutoff to the par- 
ent distribution, and the cr-dependence of the interaction 
strength (which is controlled in the PLG by the expo- 
nent a in Eq.QJ. Additionally the model, and the results 
we have provided, may prove useful as a computation- 
ally efficient test bed for anyone wishing to "tool-up" for 
simulation studies of polydispersity on other systems. 



Finally we note that experiments have been performed 
on polydisperse polymers which report that the critical 
exponents for the critical coexistence curve are Fisher- 
renormalized with respect to the standard Ising values 
[27l ]. The present model may provide a useful platform 
for elucidating these intriguing findings; this is something 
which we hope to do in future work. 
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